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In the field of ultracold atoms in optical lattices a plethora of phenomena governed by the hopping 
energy J and the interaction energy U have been studied in recent years. However, the trapping 
potential typically present in these systems sets another energy scale and the effects of the cor- 
responding time scale on the quantum dynamics have rarely been considered. Here we study the 
quantum collapse and revival of a lattice Bose-Einstein condensate (BEC) in an arbitrary spatial 
potential, focusing on the special case of harmonic confinement. Analyzing the time evolution of 
the single-particle density matrix, we show that the physics arising at the (temporally) recurrent 
quantum phase revivals is essentially captured by an effective single particle theory. This opens the 
possibility to prepare exotic non-equilibrium condensate states with a large degree of freedom by 
engineering the underlying spatial lattice potential. 

PACS numbers: 03.75.Kk, 05.30.Jp, 03.75.Hh 



I. INTRODUCTION 

The progress in the field of ultracold atoms has opened 
a new path towards directly observing quantum dynam- 
ics in the laboratory. The interference of condensates, 
Bloch oscillations, solitons and dynamics of spinor con- 
densates are examples of interesting dynamical effects, 
which, prior to their direct realization in cold atomic sys- 
tems were hard to observe and only acted as theoretical 
textbook examples [TH5]. The first step was the real- 
ization of non- or weakly interacting condensates, which 
allowed for the observation of single particle quantum 
dynamics [6-9J. However, the more interesting systems 
are those with strong inter-particle interactions, beyond 
the scope of a single particle Gross-Pitaevskii description. 
Strongly interacting systems became available using the 
capabilities of Feshbach resonances and optical lattices. 
An intriguing example for strongly correlated quantum 
dynamics is collapse and revival (CR) of a condensate 
when the lattice is suddenly ramped up and local inter- 
actions dominate [FJ [TUHT2"] . This has been successfully 
used to observe physics beyond the Bose-Hubbard model, 
quantifying the density-dependence of the intra- and in- 
terspecies interaction energy arising from the admixture 
of higher band contributions [13fH6j . In recent works 
[TTHin] the effect of a trapping potential on the n(k = 0) 
component during CR, as well as a renormalization of 
the revival time by a finite J have been considered. 

Here we study the full momentum distribution dur- 
ing CR dynamics. We show that after a sudden lattice 
quench, the dynamics in an inhomogeneous spatial po- 
tential can be described by an effective single-particle 
theory at the discrete times of quantum phase revivals 
in the limit of negligible J. For the specific case of a 
harmonic potential, condensate states consisting of co- 
herent superpositions of discrete, equidistantly spaced 
quasi-momentum states are obtained, forming an intrigu- 



ing case of non-equilibrium state preparation. The po- 
sition of the quasi-momentum components in the first 
Brillouin zone is highly sensitive to a shift of the trap- 
ping potential relative to the lattice (see Fig. [I]) and the 
appearance of the peaked momentum pattern requires 
the harmonic trapping time scale to be a rational mul- 
tiple of the revival time. In an experiment this allows 
to extract the global trapping potential with high pre- 
cision. We derive an approximate analytic theory based 
solely on the dynamical evolution of the single-particle 
density matrix (SPDM). We verify the applicability un- 
der experimental conditions by furthermore performing a 
numerical analysis based on dynamic bosonic Gutzwiller 
theory, which includes condensate depletion in the initial 
state, finite size effects and finite tunneling. 

We briefly review the basic CR physics, and proceed 
by extending the analysis to inhomogeneous systems, in- 
cluding other effects, such as finite tunneling and density- 
dependent interaction parameters. To observe CR, a 
bosonic system of atoms in an optical lattice is prepared 
in a Bose condensed state, before the lattice depth s 
(expressed in units of the recoil energy E r ) is suddenly 
ramped up [5J [T31 H] (typically to s > 25) , such that the 
interaction strength U becomes the dominating energy 
scale. 



II. HOMOGENEOUS CASE 

In the case of a homogeneous system, where all lattice 
sites are equivalent and s is sufficiently high that J/U 
can be neglected [20], the subsequent time evolution is 
generated by the Hamiltonian % = ^n(n — 1). The evo- 
lution of a local annihilation operator can be expressed 
in the Heisenberg representation as (choosing h = 1) 

b(t) = e mt be- mt = e~ lUilt b. (1) 
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III. INHOMOGENEOUS CASE 



FIG. 1. (Color online). The trap shift lo is defined as the dis- 
tance between the minimum of the harmonic potential and the 
center of the next lattice site (in units of the lattice constant 
a). As the confinement is created by laser beams different 
from those creating the lattice in a given direction, thermal 
and mechanical fluctuations in the experimental setup can 
lead to a change in 1q. 



We note that the bosonic Gutzwiller method is well jus- 
tified here, as it assumes a product state of local on-site 
states. The local number statistics in the superfluid (SF) 
before the ramp-up, including number squeezing at ini- 
tially finite U I J is well described and becomes poissonian 
in the limit U / J —> 0. In this ideal limit, the ground state 
for a system consisting of N particles on L sites becomes 
a coherent state in the k = mode \ip) = \z)i ® . . . (?) \z) l 
within a grand-canonical, [/(l)-symmetry breaking de- 
scription with a density n = N/L = \z\ 2 . For this state 
the expectation value of the order parameter (JlJ during 
CR can then be evaluated exactly [321 EH [22] 



1 p(t) = (b(t)) =zexp(|z| 2 (e- 



-iUt 
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The condensate fraction is defined as f c = Xq/N, where 
Ao is the largest eigenvalue of the SPDM p\,y = 
and 1 is the vectorial site index on the 3D lattice. In the 
limit of a large homogeneous system, Eq. ^ leads to 



= 2| 2 | 2 (cos(C/t)-l) 



(3) 



which is monotonically related to the visibility measured 
in experiments |131 123j . At times t m = m t rQV with the 
revival time t Iev — 2%/U and m G N, the coefficients 
c n {t) in the Fock representation \ijj{t))\ — Y^=o c n{t)\ n )\ 
periodically coincide after having performed mn(n— 1)/2 
rotations in the complex plane, which leads to a revival of 
the initial condensate. In the case of finite interactions 
U I J > prior to the ramp-up, the phenomenological 
ansatz 



fc(t) 



a e 



20(cos([/t)-l) 
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(4) 



leads to a remarkably good agreement with the exact nu- 
merical results for arbitrarily long times. The fitting pa- 
rameters a, 0, 7 G R incorporate number squeezing and 
depletion in the initial condensate. 



For an inhomogeneous system without translational 
symmetry, the single site description does not suffice. 
The time evolution in the atomic limit (neglecting J/U) 
is generated by the lattice Hamiltonian % = y h\{n\— 
1)— A*i"-i with the effective chemical potential [i\ at lat- 
tice site 1 accounting for the inhomogeneity. Analogously, 
the time evolution of the local annihilation operators in 
the Heisenberg picture is 



5 mt b ie - 



mt 



(5) 



In case of an initially ideal (but not necessarily homo- 
geneous) condensate, the time- and position-dependent 
local order parameter can be evaluated exactly 



(6) 



We now turn to the condensate fraction by the convenient 
decomposition of the SPDM using Eq. [5] 



Pi 



,v(t) = e-^ibl e iC/ ( b !" h i- b F b i') t 6 1 ,) e ^>' t . 



(7) 



Remarkably, the time-dependent condensate fraction 
f c (t) in an inhomogeneous system is identical to that 
in the homogeneous system (i.e. /ii = fi) for a suffi- 
ciently deep lattice (for J/U w to hold) and identical 
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FIG. 2. (Color online). Time-dependent occupation of the 
k = mode (a, b) and condensate fraction (c) for a 87 Rb 
condensate prepared in a 738nm lattice with a central density 
of n = 2.6. The lattice is ramped up from s — 8 and a 
157.4Hz trap to s — 26 and a 346.2Hz trapping frequency with 
lo,i = 0.5, such that t tr a P = I00t rC v = 28.29ms. In subplot 
(a), J is artificially set to zero after the ramp to compare the 
numerical results to Eq. [9] The inset in (b) shows the initially 
strong decline in the occupation of the k = mode (blue 
line) and the central peak N p (t) (dotted orange line), see [24] . 
Additionally, the dashed line in panel (b) shows n(k = 0, t) 
for a Gaussian trap, with parameters corresponding to the 

^ = 0.1. 



harmonic potential, and a cloud extension ratio 
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FIG. 3. (Color online). Quasi-momentum profile in the first Brillouin zone for different times and lattice shifts lo during 
CR in a harmonic trap, individually scaled on the vertical axis for visual clarity. Parameters as in Fig. [2j but initially in a 
27.5Hz trap with a central density of n — 1.67. Subsequently, this is ramped up to s = 40 and a 197.3Hz trapping frequency, 



such that ttr 



200t r , 



43.88ms. Each of the plots in the upper row depicts the projected quasi-momentum density 



n x (k x ) oc f^ a , a dk y dk z n(k) for the three 3D plots (for a given time, but different trap shifts) in the respective column. These 

projected distributions are all normalized such that dk x n(k x ) — 1, where k x — k x a/n. A trap shift of lo in the x-direction 
does not change the overall structure of n(k), but leads to a translation of the quasi-momentum profile by exactly 2W t lo in the 
A^-direction. This shift is indicated by the arrows in the projected quasi-momentum profiles in the upper row. 



U. Here, the SPDM p(t) = U(t)p(t)W(t) is connected 
to the SPDM of the homogeneous system during CR at 
the same time p hV {t) = (b\ e lUib ' b ^ b l' b ^ t b l ,) through a 
time dependent unitary transformation by a diagonal ma- 
trix with elements U\\>(t) — e -1 ^ 1 *^/ , thus leaving the 
eigenvalues and hence the condensate fraction invariant 
(we emphasize that this does not rely on any approxi- 
mation). This equivalence means that the revival of the 
condensate is actually not damped by spatial inhomo- 
geneities. At revival times the condensate state differs 
from the k = Bloch state, which appears as damping 
in certain observables such as the visibility [5l [TOj [12] , but 
is fundamentally different from a real damping of the con- 
densate. This relation also holds in the case of depleted 
condensates, where finite interactions are present prior 
to the ramp up, as also reflected by the numerical results 
for a harmonic trap shown in Fig. [2j whereas f c (t) (sub- 
plot (c)) is only damped by a finite J, n(k = 0) is almost 
completely suppressed after a single revival time by the 
very strong harmonic trap (not considering much longer 
times, where it may reappear as discussed later). 

To calculate the condensate fraction for the specific 
structure of the SPDM within a Gutzwiller state, a highly 
efficient method was found to be the iteration of the self- 
consistency condition for the largest eigenvalue Aq in the 



i-th iteration step 

A i+1) = ^ £ -TT^ 1 > (8) 

obtained from a rearrangement of the eigenvalue equation 
(see appendix A for detailed discussion). It can be shown 
that on the interval [maxj(rij — |^j| 2 ),iV tot ], the largest 
eigenvalue is the only attractive fixed point and rapid 
convergence is achieved. In the case of strong depletion, 
an adapted Lanczos algorithm is used. 

Reverting to the Schrodinger picture of the SPDM, the 
decomposed form shows that if a condensate state (i.e. a 
normalized eigenvector (/)[ ' (t) of p(t) corresponding to 
the macroscopic eigenvalue) is present at time t, its time 
evolution is exclusively determined by the operator IP if) , 

i.e. <j)f\t) = £ 1( Uf v (t)(j)!jp(0). This time evolution is 
formally also valid during collapse, where the conden- 
sate fraction is suppressed. Its observation, however, is 
only possible at revival times t m , which can be tuned 
independently, thus effectively allowing the experimen- 
tal observation of single particle dynamics in a lattice at 
arbitrary times. Note that although the effective single 
particle dynamics may resemble Gross-Pitaevskii theory, 
its derivation and validity is fundamentally different and 
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relies on the opposite limit U ^> J. 

Since the true condensate fraction is experimentally 
not directly accessible in inhomogeneous systems [251125] . 
we consider the dynamical evolution of the full quasi- 
momentum distribution n(k). Under restriction to the 
lowest lattice band, the physical momentum distribu- 
tion consists of replicas of the first Brillouin zone, scaled 
by the Wannier function's respective Fourier coefficient. 
Using the spectral decomposition of the SPDM p\,y — 

^2f = Q Xi (j)^ 4>^) in conjunction with the defining prop- 
erty of a BEC, i.e. only a single macroscopic eigenvalue 
Aq, n(k) is found to be well approximated by 



Ml 



10 



n(k,t)«^/ c (t) 



(0) 



+ C(t). (9) 



To a good approximation, the non-condensed atoms lead 
to a constant incoherent background in n(k) and f c (t) 
can be taken from Eq. [4j exploiting its independence of 
spatial inhomogeneities during CR. The background term 
C(t) — N(l — f c (t))/L ensures particle number conserva- 
tion, J2k n (k) — N. We find that all qualitative features 
arising in a full time-dependent Gutzwiller simulation are 
captured surprisingly well by Eq. [9] at all experimentally 
relevant time scales. 



IV. SPECIAL CASE: HARMONICALLY 
TRAPPED SYSTEM 




FIG. 4. (Color online). Particle number in the central quasi- 
momentum peak (see [24]) as a function of time and harmonic 
trapping strength Vt for a trap shift lo.i = 0.5. The position 
of the maxima allows for a high precision determination of Vt 
in the presence of the optical lattice. The calculation was per- 
formed for 87 Rb in a 738nm optical lattice at s = 32, includ- 
ing finite J effects and realistic density-dependent interaction 
parameters U n [13] . Insets: n(k) at points (a),(b). 



We shall now discuss the characteristics of this time 
evolution and possible applications. For the case of an 
underlying harmonic trapping potential p\ = — Vt(l — 
1 ) 2 + po, results from a numerical dynamic Gutzwiller 
calculation for a trapped condensate (diameter ss 120a 
in the x-y-plane) are shown in Fig. [2] for n(k = 0,i) 
and the full n(k, t) in Fig. [3] for certain times. The es- 
sential dynamics of n(k, t) at revivals is well captured 
within a single particle picture for the condensate state 
\tp(t)) — J2\ (although the condensate vanishes 

between revivals), where |1) is the Wannier state at site 1. 
Its time evolution is generated by U(t), since the largest 
eigenvector of p(t) is always identical at revival times. 

The condensate's contribution to n(k, t) (relevant dur- 
ing revivals) can then be expressed as 



n(k,i) 



-i' 2 ) 



with the amplitudes A i>v = <p[ 0) (Q)^ (0). The first ex- 
ponential term shows that a translation of the harmonic 
trapping potential by 1 leads to a time-dependent over- 
all translation by 2tVtlo of the quasi-momentum profile, 
which is periodic in the first Brillouin zone. The second 
term implies a temporal periodicity with a period of the 
trap time itrap = 27r/U t , up to a translation Aiml . At 
most times the last term interferes destructively, but at 



(10) 



times til™* 1 = (n + l/m)i t r ap multi-peak structures ap- 
pear in n(k) and are observable if, additionally, is a 
multiple of t Tev . In the limit of an initially homogeneous 
condensate, the discrete Fourier transform in Eq. ( 10 1 can 

be evaluated explicitly at the times t 
B) and one obtains 



n(k,t) = N f c (t) J] 



n (see Appendix 
F(q u m) 



(11) 



with F(q h m) = 6 qiinodl . (l + 5 0ym mod 2 (-l) ™ +qi ) 
and q = Tpk — 21o. This constitutes the dominant con- 



tribution in the outermost left (tg ) and right columns 

(t 2 ') of Fig. [3] during revival, whereas the peaks are 
strongly washed out in the second last column, where the 
revival time does not fully coincide with (t^). During 
collapse, shown in the second column of Fig. [3] the dis- 
tribution is essentially flat with little substructure. We 
stress that at these specific times the condensed atoms 
collectively occupy one single particle state, consisting 
of a coherent superposition of quasi-momentum states 
(K(m,is x ), K(m,Uy)) with equal weight, i.e. for the 2D 
projected case in a mode 



1 

rn 



E 



(12) 




with %K(m, Vi) = [±{2{vi + lo, t ) + 1 + f )] mod 1. 
This is not to be mistaken as a signature of a fragmented 
condensate. As the appearance of these patterns is gov- 
erned not only by the inverse trap time scale, but also by 
the CR time scale, the temporal discretization is implic- 
itly given by the revival time t rcv . A plot of n(k = 0, t) as 
a function of the trapping frequency is shown in Fig. [4] 
the 1/Vt scaling of the time at which a pattern with a 
significant n(k = 0, t) component appears is in full accor- 
dance with Eq. [lO] and the fainter hyperbolas correspond 
to momentum profiles with multiple spikes, where, cor- 
respondingly, the k — component is weaker. The non- 
persistent, dotted intensity of the hyperbolas directly re- 
flects the interplay of the two time scales: only when a 
specific trap pattern coincides with a revival of the mat- 
ter wave field, a significant n(k = 0) is observed. 

A further quantity that can be directly extracted from 
the distributions in Fig. [3] at the revival times, is the po- 
sition of the trapping potential in a scale smaller than the 
lattice spacing a. As can be seen from the term —2tV t lo 
in Eq. |10| a microscopic shift 1q in the harmonic trap- 
ping potential leads to a time-dependent translation of 
the quasi- momentum distribution n(k), as visible in the 
rows with different l 0iX = 0, 0.2, 0.5 in Fig. [3j This trans- 
lation can be experimentally observed on a macroscopic 
scale after time-of-flight expansion. Thus by comparing 
the initial ri(k, t = 0) with n(k, t m ) at revival times, lo 
can directly be determined with high precision. 

In turn, this may find application in characterizing 
the thermal or mechanical stability of the optical setup: 
as the confining harmonic potential is created by laser 
beams different from the ones generating the optical lat- 
tice in each dimension (see Fig. [TJ , thermal drifts [27] can 
lead to a shift of 1 , which can be observed in n(k, t). 



V. VALIDITY OF OUR PREDICTIONS 

A crucial question is whether our predictions are ro- 
bust enough to be observable for the harmonic trapping 
potential in typical experimental setups. As the most im- 
portant effects we consider the sensitivity towards precise 
timing as well as deviations from the idealized harmonic 
trap shape due to the Gaussian intensity profile of the 
laser beam. For a given trap strength V t , the required 
temporal precision (i.e. the time span between lattice 
ramp and the release of the cloud) is of the order of the 
revival time, as can be seen in Fig. 2 and Fig. 4. This 
time scale is precisely controlled in current experiments 
[TBI IT5] . The deviations originating from the Gaussian 
trapping potential require a somewhat more quantita- 
tive analysis. The relevant parameter, which quantifies 
the quality of the harmonic approximation is the ratio 
of the cloud radius and the beam radius. In Fig. 2(b) 
we show the time-dependent occupation n(k = 0, t) for 
a true Gaussian trapping potential for the cloud exten- 
sion ratio ~ Rc '°" d = 0.1 in comparison to the idealized 
harmonic case. We find that in this experimentally real- 




FIG. 5. (Color online). Comparison of the quasi-momentum 
profiles in the first Brillouin zone for quadratic (a) and Gaus- 
sian trapping potentials (b,c,d). The relevant parameter 
quantifying the harmonicity of the Gaussian trap is the ra- 
tio between the cloud radius and the width of the Gaussian 
potential (extension ratio). For subplots (b,c,d) this value was 
set to (0.1, 0.15, 0.2) respectively, whereas the trap strength 
and other parameters were set to agree with the strength of 
the harmonic potential in Fig. 3 in the linear approximation 
(t = 25t rcv and lattice shift lo = (0.5,0.5,0.5)). Compared 
to the harmonic trap, the peaks are successively broadened 
with an increasing filling factor for a Gaussian potential, but 
remain clearly visible, even for a large extension factor of 0.2. 
Due to the scaling of the 5-response of the peaks in n(k) for 
the harmonic case, the theoretically predicted height of such 
a peak on a single point of the discrete fc-lattice is strongly 
decreased (cf. black dotted line in Fig. 2b). In the experi- 
mental detection, this should be less visible, as the observed 
quasi-momentum distribution is in any case convolved with 
the point spread function of the optical imaging system. 



istic regime the peaks are somewhat smaller, shorter and 
broader. However, on the time scale of t trap this effect 
is weak and should not wash out the structure of n(k) 
beyond experimental resolution. This can also be seen 
in Fig. 5, where a comparison of the quasi-momentum 
profile after t = 25i rcv is given for a harmonic and Gaus- 
sian trap with different cloud extension ratios. In a very 
recent experiment [28], exactly these quasi-momentum 
patterns have been observed in the non-interacting one- 
dimensional system. In our picture, those data corre- 
spond to trap dynamics with an infinite revival time. 
Since the trapping potential in |28j is not perfectly har- 
monic, the general feasibility of our predictions has al- 
ready been experimentally confirmed. 

The creation and observation of non-equilibrium con- 
densate states in a non-interacting, one-dimensional sys- 
tem as reported in |28j is contained in our theory as a 
limiting case. The quasi-momentum distribution of a 
corresponding one-dimensional system is nothing but the 
projection of the two-dimensional distribution on the k x - 
axis. The emergent patterns in momentum space can 
be viewed as the temporal matter-wave analogon of the 
Talbot effect known from classical optics (29] [30] . From 
this point of view, one may also see our theory as a gen- 
eralization of the temporal matter-wave Talbot effect to 
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higher dimensional interacting systems (i.e. two or three 
dimensions). 



VI. CONCLUSION 

In conclusion, we have analyzed the full momentum 
distribution during collapse and revival of an atomic con- 
densate after a lattice quench, showing that the state at 
revival times can effectively be described by a single par- 
ticle theory. We have shown that the true condensate 
fraction is not damped by the presence of an arbitrary 
spatial potential, but that the condensate state takes on 
a different form, such that a damping seems to appear 
when only looking at n(k = 0,t). For the specific case 
of a harmonic trapping potential, condensate states fea- 
turing a periodic peak structure in momentum space are 
created, allowing for high precision determination of the 
total trapping frequency of the underlying optical and 
magnetic trap, as well as the position of the trapping 
potential on a sub-lattice scale. By engineering the spa- 
tial potential after the ramp-up, our concept is suitable 
for the preparation of exotic non-equilibrium condensate 
states, as well as transforming different condensate states 
into each other. These prospects seem realistic in the 
light of recent experimental progress on the engineering 
of arbitrary spatial potentials 
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Appendix A: Efficient Calculation of the Condensate 
Fraction for a Bosonic Gutzwiller State 



To obtain the condensate fraction f c in an inhomo- 
geneous system, it is necessary to determine the largest 
eigenvalue Ao of the SPDM pij> = {b\b v ), which is a her- 
mitian L x L-matrix, L being the number of lattice sites. 
The condensate fraction is then defined as f c = -^ s -, 

J JVtot 

where N tot is the total particle number. Since we deter- 
mine f c (t) at many different times, an efficient algorithm 
to determine the largest eigenvalue is desirable. Here we 
present a method to determine the largest eigenvalue of 
the SPDM of a variational state of the bosonic Gutzwiller 
form, where the computational effort scales proportion- 
ally with L. We firstly derive relations holding for all 
eigenvalues of a SPDM of this form and subsequently 
use these to formulate a fix point relation for the eigen- 
values. The largest eigenvalue can always be found by 
numerically iterating this fix point relation, as the largest 
eigenvalue on a large interval around the fix point. 



We begin by considering the explicit form of the SPDM 
for a bosonic Gutzwiller state: 



Pl , v ={b\b v ) =^*#+*i,i. (^-IV^I 2 ) 



(Al) 



the off-diagonal part factorizes into products ^ ipf of 
local order parameters tpi = (6;), whereas the diagonal 
part consists of the local densities n; = (&J&;). Now let 
v be an eigenvector of p with corresponding eigenvalue 
A (i.e. pw = Av). Focusing on a single element j of an 
eigenvector, we obtain 



A Vj = ^2 Pj,l v l = ^ i'l + ( n i 



. (A2) 

Rearranging this equation and multiplying both sides by 
tpj leads to 



I 



Finally, we divide both sides by A 
over all lattice sites j to obtain 

3 



A 



ikvt- (A3) 
f \tjjj\ 2 and sum 

(A4) 



This expression no longer contains the eigenvector v and 
therefore holds for all eigenvalues A of p. Defining a = 
maxj*{(rij — IV'jl 2 )}: one can directly infer that p cannot 
have two different eigenvalues A and A' in the interval 
(a, AT tot ] both fulfilling Eq. A4 In this interval, all the 
addends in Eq. A4 are positive and for A, A' € (a, N tot ] 
and A < A' we find that 



1 



1^ 



A — n 7 

3 J 



\4>i 



> 



^ A' — n,i 

3 J 



\i>3 



1 (A5) 



leads to a contradiction. This reveals the interest- 
ing property, that from the intrinsic form of a bosonic 
Gutzwiller state, this cannot describe a fragmented BEC. 



We now multiply Eq. A4 by A, which can then be seen 
as a fixed point, which is attractive in the interval A € 
(a,N tot ). With this equation one can define the fixed 
point iteration 



A«|^- 



At' 



(A6) 



This iteration, restricted to A 1 € (a,N tat \, always con- 
verges to the largest eigenvalue of the SPDM with ar- 
bitrary precision and directly scales with the number of 
lattice sites L. In most cases a relative error ^A < io~ 6 
is already obtained after 3-15 iterations. 



Appendix B: Explicit evaluation of the 
many-particle state at certain rational multiples of 
the harmonic oscillator time 

The aim of this section is to evaluate |F(n)| 2 , physi- 
cally corresponding to the quasi-momentum distribution 
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of the condensate state at certain points in time after 
evolving in a harmonic confining potential, with 



L-l 



F(n) = ]T 



(Bl) 



1=0 



Here n, L £ N are integers and L is an integer multiple 
of to, i.e. L/m £ N, to avoid broadening of the peaks in 
the resulting spectrum. 

The summation integer I is split into I = p + qm, 
where the newly defined integers can independently take 
on values in the ranges p £ {0, ...,m — 1} and q £ 
{0, — 1}. Thus, the previous sum can be written 
as 



m-i sr-i 

p=0 q=0 

= £ 

Using the identity 

AT-l 



27rif (p+qrn) (p+ ^" ) 



m — 1 

£- 

p=0 



(B2) 



n=0 



(B3) 



for AT — — . one obtains 

m ' 

L „ 



F(n) 



n mod 



i 53 «KW 

p=0 



(B4) 



This restricts the resulting spectrum to a form with 
non-zero values only at to periodically spaced positions 
n = r m w ith r £ N. To determine the value of the sum 
at these discrete values of n, it is convenient to define the 
function 



9( r ) = £ £ 

p=0 



(B5) 



and it will now be shown that | <? C 7 ™ ) 1 2 * s independent of 
r <E N for a fixed value of m. 



rn— 1 

lsMI 2 = £ e 

p,p'=0 

m — 1 p— 1 

= m + 

p=l p'=0 



Tl— 1 p— 1 



(B6) 



We now perform a transformation of the summation in- 
dices and define d = p—p' and s = p+p'. The summation 
over all terms can be reexpressed in the new variables as 



m— 1 p— 1 m— 1 2(m— 1) — d 

££-£ £ ■ 



(B7) 




FIG. 6. Illustration of the resummation procedure. The 
sum of the diagonal values amounts to m, while the values in 
the lower right triangle correspond to the complex conjugated 
values in the upper left triangle, mirrored on the diagonal. A 
summation over s = p + p' at fixed d = p — p' can be under- 
stood as sums of the dth lower off-diagonal, where subsequent 
values of s are spaced by 2. 



Subsequently, 



m-l 2(m-l)-d 

| 5 (r)| 2 =m+53 53 [e 2 -^ 



m — 1 



(B8) 



53 [e 2 ™^ h(d,m)+h.< 



d=l 



where we defined the function 



h(d, to) — 5^ e 



2ni- 



(B9) 



with the summation variable q = s/2 now running over 
values with a step size of 1. To evaluate this function for 
fixed to, the explicitly known expression for the geometric 
series can be used as long as d ^ ^. I n the latter case, 
it can easily be evaluated and Eq. |B9| takes on a value 
h(d = m/2, m) = (— l) d d, leading to the generally valid 
expression 



h(d,m) = (1 - 5 d;m / 2 ) 

= (1 - bd,m/2) 



3 47ri^(m-|) 1 _ f , 



4TTI A # 



- +<5 d , m /2(- 1 ) d rf 

1 — e m 1 — K m / 

'2e- 2?rl s sin(27r^)sin(27r^)\ . 



cos(4^) - 1 



+ <5d,m/a(-l) d. 



(BIO) 



r 



Inserting this result into Eq. B8 one obtains 



\g{r)\ 2 = m + <5 , m mod 2 (-1) * m cos(7rr) + ^ (1 - S djm / 2 ) 



'4cos(27r^(r- 1)) sin(27r^) sin(27r£)' 
. cos(4^)_r ~ 



(Bll) 



The second term only contributes for even m, whereas 
the sum is identically equal to zero, as will now be shown. 
The term (1 — <5d,m/2) ensures that only an even number 
of terms may contribute in the sum. This allows for a 
regrouping of terms and a corresponding resummation, 
where the terms d = c and d — m — c are considered 
as pairs, subsequently amounting to zero by using the 
periodicity and symmetry of the trigonometric functions: 



\g(r)\ 2 =m + 5 , m mod 2 (-1) 



l -+r 



c=l 



4cos(27r^(r- 1)) sm(27rM sin(27r^) 

C0S(47T — ) - 1 

4cos(27r^(r- 1)) sin(27r^^) sin(27r^)\ 
+ cos(4tt^) ^1 J 

= m(l + <5 . m mod 2 (-l)T+' r ). 

(B12) 

We draw the connection to the quasi-momentum dis- 
tribution from the condensate state in a harmonic trap 
at time t 

n ( k)t ) = 1^^(0)^,(0)* e -*(»*-2tv t j )(J-0 e ^tv ti i 2 -i' 2 



LI' 



VL l 



)i e -itv t v 



(B13) 



where the initial condensate state is (f>i(0). Assuming 
that the condensate is initially in a homogeneous state, 
i.e. <pi(0) = 0(0) = e lVo L~z with an irrelevant overall 
phase <po, the comparison of Eq. B13 with F(n) shows 



that at the specific times t = by equating 

4ir n 

ak to = 2ir — . 

m L 



(B14) 



and using r = ™ , one can relate 



n(k,t) 



L 

\m\ 2 l 2 



L m 



F In 

6 



akl 211q 
2tt m 

' akm 



2 n mod — ,0 
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2n 



2Zr 



Thus, at the specific time t = 



X (l+<5 ,mmod2(-l)^ + ^" 2i °) 



(B15) 



(B16) 



If the condensate is initially not in a perfectly homo- 
geneous state, but varies slowly in density, the peaks are 
broadened, but the underlying structure still remains. 

For odd integer values of m, the term in brackets al- 
ways takes on the value one and the momentum pro- 
file consists of an equidistant series of peaks at k = 
— (r + 2l ) and r£Z. As required from the unitarity 
of the Fourier transform, the quasi-momentum profile is 
of course normalized ^ fe n(k, t — = 1 a t these cer- 
tain times and each peak carries exactly the same weight 
1/m. 

For even integer values of m, Eq. |B16| can be written 



as 



n(k ' i) = V™' 5 (T[l+^]-'o)modl,0 



(B17) 



and the quasi-momentum profile consists of only m/2 
equally strong and equidistantly spaced peaks at k — 

1 t^( 5 + 'o)-i] withsez. 

This analytical result agrees perfectly with the results 
of our numerical simulations. To the best of our knowl- 
edge, this analytical expression, which is also of direct 
relevance for the Talbot effect in quantum optics, has 
not been obtained so far. 
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